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Abstract 

This paper presents an investigation of the dynamics of two coupled non-identical FitzHugh-Nagumo neurons with 
delayed synaptic connection. We consider coupling strength and time delay as bifurcation parameters, and try to 
classify all possible dynamics which is fairly rich. The neural system exhibits a unique rest point or three ones 
for the different values of coupling strength by employing the pitchfork bifurcation of non-trivial rest point. The 
asymptotic stability and possible Hopf bifurcations of the trivial rest point are studied by analyzing the corresponding 
characteristic equation. Homoclinic, fold, and pitchfork bifurcations of limit cycles are found. The delay-dependent 
stability regions are illustrated in the parameter plane, through which the double-Hopf bifurcation points can be 
obtained from the intersection points of two branches of Hopf bifurcation. The dynamical behavior of the system may 
exhibit one, two, or three different periodic solutions due to pitchfork cycle and torus bifurcations (Neimark-S acker 
bifurcation in the Poincare map of a limit cycle), of which detection was impossible without exact and systematic 
dynamical study. In addition, Hopf, double-Hopf, and torus bifurcations of the non trivial rest points are found. 
Bifurcation diagrams are obtained numerically or analytically from the mathematical model and the parameter regions 
of different behaviors are clarified. 

Keywords: FitzHugh-Nagumo neural model. Delay differential equation, Double-Hopf bifurcation, Hopf-Pitchfork 
bifurcation. Torus bifurcation. Fold of limit cycles 
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1. Introduction 

Neurons and their interactions are generally assumed 
to be the determinant of the brain performance. The 
simplest model to display features of neural interac¬ 
tions consists of two coupled neurons or neural systems. 
Starting from such simple and reduced networks, larger 
networks can be built and their features may be stud¬ 
ied. In order to study complicated interaction between 
neurons in large neural networks, the neurons are of¬ 
ten put up into highly connected sub-networks or syn¬ 
chronized sub-ensembles. Such neural populations are 
usually spatially localized fTl. In this way, the model 
of two mutually coupled neurons may also apply as a 
framework of two coupled neural sub-ensembles. 
Destexhe et al. 0, by investigating two interconnected 
thalamic spindle oscillatory neurons, showed that how 
more complex dynamics emanates in ring networks with 
nearest neighbors and fully reciprocal connectivity, or 


*CoiTesponding author: Tel.:+982166165626. 
E-mail address: farajzadeh@mehr.sharif.ir 


in networks in which every neuron connects to all other 
nearby neurons. Zhou et al Q modeled a neural net¬ 
work as a small sub-network of excitable elements to 
study synchronization dynamics and the hierarchically 
clustered organization of excitable neurons in a network 
of networks. 

The FitzHugh-Nagumo (FHN) model with cubic non¬ 
linearity which was derived as a simplified model of the 
famous Hodgkin-Huxley (HH) model ID by FitzHugh 
m and Nagumo et al. 0, is a classic oscillator ex¬ 
hibiting variety of nonlinear phenomena in planar au¬ 
tonomous systems. The FHN-like systems are of fun¬ 
damental importance for describing the qualitative na¬ 
ture of nerve impulse propagation and neural activity. 
In fact, this model seems reach enough, and can cap¬ 
ture neural excitability of original HH equation Q. In 
this way, Guckenheimer and Kuehn 0 investigated ho¬ 
moclinic orbits of the FitzHughNagumo equation from 
the view-point of fast-slow dynamical systems. Hoff 
et al. ||9l investigated unidirectional and bidirectional 
electrical couplings of two identical FHN neurons, in¬ 
volving the effects of the four parameters on bifurca- 
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tions and synchronization. In modeling studies of trans¬ 
membrane potential oscillations, Bachelet et al. cni ex¬ 
plained that bursting oscillations appear quite naturally 
for two coupled FHN equations. Their coupling model 
is quite different from other investigations. Different fir¬ 
ing patterns such as chaotic firing is found in a pair of 
identical FHN elements with phase-repulsive coupling 

Ca¬ 
in addition to study of coupled neurons, the dynamics 
of a network of FHN neurons is also studied by re¬ 
searchers. Oscillatory chaotic behaviors and the mech¬ 
anism of dynamic change from chaotically spiking to 
firing death is induced in ifT^ . by considering the cou¬ 
pling of FHN oscillators in a network. Cattani C3 
proposed a network model of instantaneous propagation 
of signals between neurons in which the neurons need 
not be physically close to each other. The dynamics of 
the action potential and the recovery variable within a 
neural network is described by proposed network. The 
oscillation patterns in a symmetric network of modi¬ 
fied FHN neurons is discussed by Murza lfT4l . The au¬ 
thor described the building block structure of a three- 
dimensional lattice shaped as a torus, and identified the 
symmetry group acting on the electrically coupled dif¬ 
ferential systems located at the nodes of the lattice. 

It is known that signal transmission in coupled neurons 
is not instantaneous in general CSl. Hence a time de¬ 
lay can occur in the coupling between neurons or in a 
self-feedback loop. For example, synaptic communi¬ 
cation between neurons depends on the propagation of 
action potentials across axons. The finite conduction 
speed and the information processing time in synapses 
lead to a conduction delay. The speed of signal trans¬ 
mission through unmyelinated axonal fibers is in the 
order of 1 m/s, which results in existence of time de¬ 
lays up to 80 ms for propagation through the cortical 
network Ill5l[l6l. Moreover, the interaction delay be¬ 
tween the oscillators can be as large as the oscillation 
period ini. Therefore, time delay is inevitable in sig¬ 
nal transmission for real neurons. On the other hand, 
time-delayed feedback mechanisms might be intention¬ 
ally implemented to control neural disturbances, e.g. 
to suppress alpha rhytm IfTSl and undesired synchrony 
of firing neurons in Parkinson’s disease or epilepsy 
HSlIlOlEI]. Various delayed feedback loops have been 
proposed as effective and powerful therapy of neuro¬ 
logical diseases which are caused by synchronization of 
neurons 112011211 l22l l23l 1^ l25l . 

From the modeling viewpoint, delay in dynamical sys¬ 
tems is exhibited whenever the system’s behavior is de¬ 
pendent at least in part on its history. There are many 
technological and biological systems which exhibit such 


behavior. Some examples of delayed systems are cou¬ 
pled lasers ll2^ . high-speed milling, population dynam¬ 
ics, epidemiology im, and gene expression, 

Krupa and Touboul ||28l investigated a self-coupled de¬ 
layed FHN system to understand the effect of introduc¬ 
ing delays in a system whose dynamics, in the absence 
of delays, is characterized by a canard explosion and re¬ 
laxation oscillations. They have shown that delays sig¬ 
nificantly enrich the dynamics, leading to mixed-mode 
oscillations, bursting, canard and chaos. 

The research for coupled FHN systems with delay has 
attracted many authors’ attention. Dahlem et al. ll^ 
studied an electrically delay-coupled neural system. 
They showed that, for sufficiently large delay and cou¬ 
pling strength bi-stability of a fixed point and limit cy¬ 
cle oscillations occur, even though the single excitable 
elements display only a stable fixed point. Buric and 
Todorovic m, and also Jia et al. EH investigated 
Hopf bifurcation of coupled FHN neurons with delayed 
coupling, they concluded that the time delay can in¬ 
duce Hopf bifurcation. Buric et al. Ga, observed dif¬ 
ferent synchronization states in a coupled FHN system 
with delay and by variation of the coupling strength and 
time delay. Also they showed that the stability and the 
patterns of exactly synchronous oscillations depend on 
the type of excitability and type of coupling. In ll^ 
by employing the center manifold theory and normal 
form method, the Fold-Hopf bifurcation in a coupled 
identical FHN neural system with delay is investigated. 
The synchronization transitions between inphase and 
antiphase oscillations for delay-coupled FHN system is 
studied in If34l . Wang et al. llTSl reported that time 
delay can control the occurrence of some bifurcations 
in two synaptically coupled FHN neurons. However, 
only codimension-1 bifurcations are discussed in these 
papers. Fan and Hong 13^ . and also Xu et al. lElI, 
considered the stability and Hopf bifurcation of dou¬ 
ble delay coupled FHN system with different coupling 
strength. Yao and Tu |[38l discussed the combined effect 
of coupling strength and multiple delays on the stabil¬ 
ity of the rest point, and they obtained stability switches 
in the coupled FHN neural system with multiple delays. 
In addition to the study of Hopf bifurcation, Zhen and 
Xu 133, and also Rankovic ioi discussed the effects 
of the coupling strength and delay on steady state bifur¬ 
cations in coupled systems with identical neurons. The 
Bautin bifurcation due to small delay, was analyzed by 
Buric and Todorovic ED. Nevertheless, in that paper 
the time-delayed argument of the coupling function was 
replaced by its Taylor expansion, so that their delay- 
differential system was approximated by ordinary dif¬ 
ferential equations. Zhen and Xu Il42l studied Bautin bi- 
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furcation of completely synchronous FHN neurons. Li 
and Jiang ll43l . investigated Hopf and Bogdanov-Takens 
bifurcations in a coupled FHN neural system with gap 
junction for identical neurons. In Il44ll it was shown that 
in a simple motif of two coupled FHN systems with two 
different delay times and with self-feedback, complex 
dynamics arise. In particular, resonance effects between 
the different delay times proved to be crucial. 

According to above discussion, FitzHugh-Nagumo 
model with delayed coupling is interesting and exhibits 
rich dynamical behavior in comparison with a simplified 
version of FitzHugh-Nagumo equation ll45l . Indeed, 
the dynamics of a pair of excitable systems with time- 
delayed coupling is quite different from the dynamics 
with the instantaneous coupling. It is also known that, 
the behavior of coupled delayed systems differs from 
that of coupled oscillators with the same coupling ll35]l . 
In addition to study of coupled FHN neurons, the dy¬ 
namics of a generalized model consisting of N sets 
of FHN equations has been investigated by Buric and 
Todorovic If30l . In their analysis, they applied the ter¬ 
minology found in the description of a small lattice net¬ 
work with a finite time delay between sections for trans¬ 
ferring an excitation signal from one element to another. 
For the case of A = 2, they studied the properties of 
stability and the bifurcation of the signals. Kantner et 
al. Il46l showed that how a variety of stable spatio- 
temporal periodic patterns can be created in 2D-lattices 
of coupled oscillators with non-homogeneous coupling 
delays. The effects of heterogeneous coupling delays in 
complex networks of excitable elements described by 
the FHN model is studied in PTi . The effects of dis¬ 
crete as well as of uni- and bi-modal continuous distri¬ 
butions are studied with a focus on regular, small-world, 
and random networks. The appearance and stability of 
spatio-temporal periodic patterns in unidirectional rings 
of coupled FHN systems interacting via excitatory de¬ 
layed chemical synapses is described in Il48]l . A network 
of connected FHN neurons with delayed coupling and 
different synaptic strength of self-connection is studied 
in 1491 . The stability and oscillation of the solutions is 
investigated by constructing of Liapunov functional and 
applying the Chafees limit cycle criterion. Zhen and 
Xu and also Lin ll50l generalized the FHN models 
to a three coupled FHN neurons and Feng and Plam- 
ondon ED to a four coupled FHN neurons with time 
delay. They derived some criteria to determine the pe¬ 
riodic oscillation were provided in the multiple delayed 
FHN neural system with nonidentical cells. 

In this paper, we investigate the effect of the coupling 
strength and time delay on the stability and bifurcations 
of the system of two coupled FHN neurons. Note that 


neurons in the work of Zhen and Xu ll39l . and also Li 
and Jiang ll43l are assumed identical which leads to a 
symmetry for their system of equations, that can not be 
found in our system. Instead the symmetry in our solu¬ 
tions is due to the inherent symmetry of the FHN neu¬ 
rons, namely our system of equations is equivariant un¬ 
der antipodal map. As a result of the inherent symmetry 
of the FHN model some pitchfork, pitchfork cycle and 
Hopf-pitchfork bifurcations are observed in our model. 
We can see that the time delay can either suppress peri¬ 
odic spiking or induce new periodic spiking, depending 
on the value of the time delay. While in many of the 
above mentioned papers, only local and codimension- 
1 bifurcations of FHN systems are studied, the study of 
global and codimension-2 bifurcations are crucial to un¬ 
derstand the mechanism of different dynamical behav¬ 
iors such as transition between synchronous and anti¬ 
phase solutions and simultaneous existence of stable pe¬ 
riodic solutions. Actually this kind of bifurcation study 
is a serious task for the prediction and detection of dif¬ 
ferent dynamical behaviors and special types of solu¬ 
tions. For instance detecting multiple stable periodic 
solutions which exist simultaneously is rare in litera¬ 
ture. In this work we couldn’t detect the third stable 
periodic solution just with simulations. Actually, by the 
guidance of bifurcation study, we first predicted it’s ex¬ 
istence, then approved it with simulation. 

To study the delay effects on neural system in details, in 
this paper, we have analyzed double-Hopf bifurcations 
of our system, which is a codimension-two singularity. 
Various dynamical behaviors are classified in the neigh¬ 
borhood of double-Hopf points. We also show the exis¬ 
tence of one, two, or three different modes of spiking as 
a result of pitchfork cycle and torus (Neimark-Sacker) 
bifurcations. We want to emphasize that, by Bendixson- 
Dulac theorem, it can be seen that the single neuron of 
the FHN model does not admit periodic solutions for 
the ranges chosen in our paper, unlike the system in Il28l 
whose dynamics, in the absence of delays, is character¬ 
ized by a canard explosion and relaxation oscillations. 
Thus the periodic solutions are due to coupling of the 
neurons, or time delay. Synchronized solutions in many 
works in the literature are caused by the symmetry of 
the proposed systems, like the study of Zhen and Xu 
and also Li and Jiang ll43l . In this paper, we show 
the existence of both almost synchronized and almost 
anti-phase solutions for non-identical neurons in differ¬ 
ent range of parameters. We will show that almost syn¬ 
chronized and almost anti-phase solutions alternatively 
appear by changing the time delay. 

The paper is organized as follows. A brief description of 
the FHN neuron model is introduced in Section |D Also 
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the stability of the trivial rest point is investigated, and 
the critical values of Hopf and pitchfork bifurcations are 
driven. In Section with the aid of numerical sim¬ 
ulations, we have obtained bifurcations for the system 
with and without delay. Finally, in Section [^numerical 
simulations are carried out for determining the bifurca¬ 
tions for a large period of time delay, and delay driven 
dynamics are analyzed. A conclusion is presented in 
Section 0 

2. Model description and basic results 

In order to consider the effects of delay in the signal 
transmission between the neurons, we use the paradig¬ 
matic FHN model, 0161 with time delay, investigated 
by Wang et al. 1351. We here use a coupled nonidenti¬ 
cal FHN neural system at which all the parameters are 
assumed to be non-negative; 

vi(f) = -V[(0 -H avi(f) - wi(f) -I- ctanh(v 2 (f - t)), 

= vi(f) - hiWl(f), 

viit) - + av 2 (t) - W 2 (f) -I- ctanh(vi(f - t)), 

W2(t) = V2(t) - b2W2(t), (1) 

where a, bi and b 2 are positive constants, vi represents 
the membrane potential, wi is a recovery variable, c 
measures the coupling strength, and t > 0 represents 
the time delay in signal transmission. We suppose that 
the function which describes the influence of a neuron 
on the other one, at time f, depends on the state of the 
neuron at some earlier time f - r. We consider c and r 
as bifurcation parameters. 

2.7. Linear stability 

The equations defining the system 0 have Z 2 - 
symmetry (invariance under the change of the sign in 
the state variables). Due to the Z 2 -symmetry, the origin 
is always a rest point of the system Q, which we call 
the trivial rest point. The characteristic equation corre¬ 
sponding to linearization of the system 0 at the trivial 
rest point is 

P(c, t) = + AA^ +BA^ + CA + D (2) 

-E{A + by){A + b2)e-^^^ 

where A = b\+b2-2a, B = bib2-2a(bi+b2)+a^+2, C - 
{a^ + \){b\+b2)—2ab\b2—2a,D — a^bib2—a(bi+b2)+l, 
and E - c^. 

Due to the presence of the delay, Eq. (|^ has infinitely 
many solutions; however, the stability of the fixed points 
is determined by a finite number of critical roots with 
largest real parts. 


2.2. Elopf bifurcation 

We want to obtain some conditions to ensure that the 
system ([T]) undergoes a single Hopf bifurcation at the 
trivial rest point (0,0,0,0), when t passes through cer¬ 
tain critical values. Substituting A = io) into Eq. Q and 
separating the real and imaginary parts, we obtain 

( i / - Ba / + D - Ei - a / cos(2tw) + b \ b2 cos(2tw) 
-I- (7>i + 7>2)ws/n(2Tw)) = 0, 

-A(x? + Ccj - E((a/ sin(2Tw) - bib 2 sin(2Tw) 

-I- (7>i + b2)cocos(2Ta))) - 0. (3) 

Eliminating t from Eq. 0 gives 

w* -r Pw® 4- QOJ^ +Rto'^ + S ^0, (4) 

where P = -2B + A^,Q^B^ + 2D- 2AC - E^, R ^ 
-2BD + C^- E\b] + bl), and S ^ - E\bib 2 f. Eor 

z = we get 

+ Pz^ + Qz^+Rz + S =0. (5) 

Since the form of Eq. 0 is identical to that of Eq. (2.4) 
in the paper of Li and Wei ll52l . we may apply Lemmas 
2.2 and 2.3 analogously. The proofs are omitted. 

Lemma 2.1. If S <0, then Eq. 0 has at least one 
positive root. 

Denote h{z) - z^ + Pz^ + Qz^ + Rz + S, then we have 
h'(z) = 4z^+3Pz^+2Qz+R. Sst4z^+3Pz^+2Qz+R = 0. 
Let y - z+^, then the Eq. ^ becomes +P\y + Q\ - 
0, where Pj = 2 _ and Gi = ^ ^ + f ■ Dehne 


A = (f 


-l+iV3 

2 ’ 


yi 

= f¥4 

Va-h^ 

-f-VA, 

y2 

= f¥4 

VAeH--( 

H - 

y3 

= f¥4 

VAe^ - 1 - 

f f - Va, 

Let Zi - 

yi-j, (*=1,2, 

,3). 



Lemma 2.2. Suppose that S > 0, then we have the 
following results. 

(i) If A > 0, then Eq. 0 has positive roots if and only if 
Zi > 0 and h{z\) < 0. 

(ii) If A < 0, then Eq. 0 has positive roots if and only if 
there exists at least one z* e {zi,Z 2 ,Z 3 } such that z* > 0 
and h(z*) < 0. 

Suppose that Eq. Q has positive roots. Without loss 
of generality, we assume that it has four positive roots, 
denoted by Zj(k = 1,2,3,4). Then Eq. Q has four 


4 









positive roots, say w, = ^/i*, i - 1,2,3,4. By Eq. (j^ 
we have 

. - Bail + D)(bi + b2)a>]^ 

sin(2TtL)) = -;—r---— 

E[{bi + b2)^a)l + {ail ~ ^ 1 ^ 2 ) 1 

(AojI - Ca)k){-ail + bib2) 

E[{bi + b2)^ajl + (0)1 - bib2)^ ’ 

X {(^t - Baal + D)(bib2 - 0)1 

COS(2t6j) = - ;—^^-— 

E[{bi + b2fa)l + {ail - ^ 1 ^ 2 )^] 

{~^^l + Ca)k){{bi + b2)a)k) 

*E[{bI + b2)^a)l + {0)1 - 

Thus, denoting 

, {a)l — Bail "*■ B)){bi + b2)a)k 

E[{bi + b2fa)l + {a)l - bib2f{ 

{Aa)l - Ca)k){-a)l + bib2) 

E[{bi + b2)^a)l + {ail ~ ^ 1 ^ 2 )^] ’ 

+ D){bib2 - 0)1 

~ E[{bi + b2fa)l + {0)1 - bib2f] 

{-Ao)l + Ca)k)({bi + b2)a)k) 

E[{bi + b2fail + {0)1 - bib2f} ’ 


Aj) 


2 ^(arccos b* +2 jn), a* Js 0, 

2 ^ (27r - arccos b* + 2 jn), a" < 0, 
{k= 1,2,3,4,; = 0,1,2,...) 


( 8 ) 


then ±ia>k is a pair of purely imaginary roots of Q with 


Jj) 

‘k ■ 


Applying Lemmas 2.1 2.2 the Routh-Hurwitz (R-H) 
criterion, and Ruan and Wei’s result lf53l according to 
Rouches theorem, we have the following results. 


Lemma 2.3. Assume that A > 0, A{B-E) > C-E{bi + 
b 2 ), D > Ebib 2 , and[C-E{bi+b 2 )][A{B-E)-C+E{bi + 
^ 2 )] > A^{D - Ebib 2 ) are satisfied, ((R-H) hypothesis). 

(i) If one of the followings holds: (a) S < 0; (b) S > 
0,D > 0,zi > 0, and h{z\) < 0; (c) S > 0, D < 0, 
and there exists z* e {zi,Z 2 ,Z 3 } such that z* > 0 and 
h{z*) < 0, then all roots o/Q have negative real parts 
when T e [ 0 , tq), such that tq = min{Tf\Tf\Tf\ ). 

(ii) If the conditions (a)-(c) of (i) are not satisfied, then 
all roots o/Q have negative real parts for all t > 0. 

Motivated by Lemma 2.3 of the work of Ruan and 
Wei II 53 I . Lemma 2.4 of the work of Li and Wei ll52]| . 
Lemma 2.5 of the paper of Hu and Huang ll54l . and 
Theorem 2.1 of Lan and Hong 13^ . we obtain following 
conclusions. 


Lemma 2.4. Suppose h'{zo) 0. Ifr — TQ, then ±ia)o 
is a pair of simple purely roots of Eq. In addition, 
0, and the sign of is consis¬ 

tent with that ofhfzf). 


Applying Lemmas |2.3|2.4 we obtain the following 
theorem immediately. 


Theorem 2.1. Suppose hypothesis (R-H) of Lemma \2.3\ 
hold. 

(i) If non of the conditions (a) S < 0; (b) S > 0,D > 
0,zi > 0, and h{z\) < 0; (c) S > 0, D < 0, and there 
exists a z* & {zi,Z 2 ,Z 3 } such that z* > 0 and h(z*) < 0 is 
satisfied, then the zero solution of0 is asymptotically 
stable for all t > 0. 

(ii) If one of the conditions (a), (b), or (c) of (i) is satis¬ 
fied, then the zero solution of system Q,/orT e [0,To), 
is asymptotically stable ( tq is the parameter defined by 
the Lemma |Z^. 

(Hi) If one of the conditions (a), (b), and (c) of (i) 
is satisfied, and h'{zl) ^ 0, then for r = ~ 

1,2,3,...), the system 0 undergoes a Hopf bifurcation 
at (0,0,0,0). 


2.3. Pitchfork bifurcation 


Now we want to study the possible steady state bifur¬ 
cations of the trivial rest point. Motivated by Lemma 2.1 
of the work of Li and Jiang l43l . we have the following 
conclusions on the eigenvalues of the system 0- 



2ahib2+2a-{a^+l-c^){bi+b2) _ c^-C 


2 c^h\b 2 2c~b\ bo ' 

A — 0 is a double root of 0, if and only if c^ — 

a^b\b2+\-a{{b\+b2)) _ D _ 2ab\b2+2a-{a^-\-\-c^){b\f-b2) _ 

b\b2 ~ b\b2’ ~ 2c^b\b2 ~ 


c-~C 

2c^b\b2 


and 


^*-fr^U=o = ^1^2(1 - 2 PP) + + 2 + (foi + b2){2TC^ - 2 a). 


We can see that our model 0 have Z 2 -symmetry, and 
in the presence of this symmetry the pitchfork bifur¬ 
cation becomes generic isn. Therefore, possible zero 
eigenvalue of the Eq. 0 is correspondent to pitchfork 
bifurcation, and saddle-node bifurcation is impossible. 
Moreover, possible double zero eigenvalue of the Eq. 
0 is correspondent to Hopf-pitchfork bifurcation. 

In the rest of the paper we will explain coupling and de¬ 
lay driven dynamics of the system 0. Bifurcation dia¬ 
grams are obtained numerically or analytically from the 
mathematical model and the parameter regions of differ¬ 
ent behavior are clarified. Numerical simulations using 
the bifurcation analysis software DDE-Biftool 15^ are 
carried out to illustrate the validity of the main results. 
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3. The dynamics for small time delays 


Bifurcation diagram (Ampiitude of vl{t) and vl of rest point) 


We want to consider possible rest points and limit cy¬ 
cles of the system ([T]), and their bifurcations for the vari¬ 
able parameters c and small t’s. As stated in the previ¬ 
ous section the origin is always a rest point of the sys¬ 
tem, without any constraint on the parameters. In the pa¬ 
per we choose the value of the parameters as a = 0.55, 
bi = 1.128, and b 2 = 0.58. This choice of values makes 
the trivial rest point stable for c = 0. Hence, by this 
choice of the parameters the neurons are at rest, without 
coupling. First, we want to study the coupling strength- 
driven dynamics of the system. Therefore, as the next 
step, let us briefly consider the system Q, for T = 0. 
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3.1. Instantaneous coupling 

We know that (0,0,0,0) is always a rest point of the 
system and according to the Z 2 -symmetry of the 
system, if (Vj, Wj, v^, wp is a rest point of the system 
Q, so is -(Vj, Wj, vp wp. Therefore, pitchfork bifur¬ 
cation is generic for our model, and the double zero 
eigenvalue of the Eq. Q is correspondent to pitchfork 
bifurcation. For systemTQ we can see that for fixed pa¬ 
rameters a - 0.55, b\ - 1.128, b 2 - 0.58, and for the 
parameter c - 0.6285, the pitchfork bifurcation occurs 
and a pair of new rest points appear. To fully study the 
bifurcations of the trivial rest point we will increase the 
parameter c continuously. Fig. [T] If we consider two 
coupled neurons without delay, for small amount of the 
parameter c the trivial rest point is stable. By increasing 
the parameter c, it is seen that for c - 0.3974 a super¬ 
critical Hopf bifurcation occurs. Fig. [^sup-H), and the 
neuron’s stability switches from resting state to periodic 
spiking. The amplitude of the stable limit cycle grows 
continuously with c, in a small interval of the values of 
c. By increasing the parameter c, there is a sub-critical 
pitchfork bifurcation at the parameter c - 0.6285, Fig. 
[3P)> and a pair of antipodal rest points appear. The 
stability of the origin before and after the pitchfork bi¬ 
furcation and the two new rest points depends on the 
occurrence of Hopf bifurcation for the trivial rest point. 
Since the super-critical Hopf bifurcation occurs before 
pitchfork bifurcation, all the three rest points of the sys¬ 
tem are unstable after pitchfork bifurcation, and the sys¬ 
tem is mono-stable, i.e. the only stable state of the sys¬ 
tem is the limit cycle that stated above, see Fig. 0 By 
increasing the parameter c, for c - 0.9751, for two non¬ 
trivial rest points a sub-critical Hopf bifurcation occurs, 
and two small unstable limit cycles appear around non 
trivial rest points. Fig. [iJsub-H). Hence, two nontriv¬ 
ial rest points become stable and the system becomes 


Figure 1: The squared (black) line shows the stable rest points of cou¬ 
pled neurons without delay. For c < 0.3974 the origin is the only 
rest point. The cross (blue) line shows the amplitude of the periodic 
cycle which appears after super-critical Hopf bifurcation (sup-H) of 
the trivial rest point. After Hopf bifurcation the trivial rest point be¬ 
comes unstable, which is shown by triangle marker. When pitchfork 
bifurcation (P) occurs at c = 0.6285, a pair of nontrivial rest points 
appear which are unstable and are shown by triangle marker. The 
nontrivial rest points become stable due to sub-critical Hopf bifurca¬ 
tion (sub-Hopf) at c = 0.9751, and two unstable limit cycles appear. 
Two unstable limit cycles collide at c = 1.0545, and homoclinic bifur¬ 
cation occurs (h). The stable and unstable limit cycles annihilate each 
other and a fold bifurcation of limit cycles (fc) occurs at c = 1.0721. 


dichotomous; two rest points and a limit cycle. By in¬ 
creasing the parameter c, two unstable limit cycles grow 
continuously until they collide at c = 1.0545, Fig. [^h). 
Therefore, as a global bifurcation, an eight-figure sad¬ 
dle homoclinic bifurcation occurs. At the bifurcation 
point the period of the unstable limit cycles have grown 
to infinity, see Fig. |^Feft). After this bifurcation a big 
unstable limit cycle around the rest points appears. We 
should emphasize that after this bifurcation the system 
remains dichotomous and the behavior of the neurons 
doesn’t change. It’s not hard to check that the saddle 
quantity is positive and according to the theorems on 
homoclinic orbits developed by Shilnikov ll57l . the bi¬ 
furcation results in the birth of a unique unstable limit 
cycle from homoclinic orbit. The unstable eight-figure 
saddle homoclinic orbit is shown in Fig. |^Right). By 
a small change of the parameter c, at c ^ 1.0721 a 
fold of limit cycles bifurcation occurs, point (fc) in Fig. 

Therefore, the stable limit cycle emanating from the 
Hopf bifurcation of the trivial rest point and the new 
born unstable limit cycle approach and annihilate each 
other, see Fig. Thus the only stable states of the sys¬ 
tem are non trivial rest points and the system becomes 
quiescent. 
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Figure 2: (Left) The period of the unstable limit cycle which appears 
through the sub-critical Hopf bifurcation. At the bifurcation point 
c = 1.0545, the period of the unstable limit cycle has grown to infinity. 
(Right) The unstable eight-figure saddle homoclinic orbit. 




t t 


Figure 4: The impact of the time delay on the dynamics of the neurons 
for c = 0.5. (a) T = 0.1, the neurons spike with a large amplitude, (b) 
T = 0.25, the neurons spike with a small amplitude and are near the 
super-critical Hopf bifurcation, (c) t = 0.4, the neurons are quies¬ 
cence. (d) r = 1.2, the neurons ai‘e spiking. 


Figure 3: Stable (blue solid line) and unstable (red dashed) limit cy¬ 
cles near the fold of limit cycles bifurcation. 


3.2. Delayed coupling 

Now we come back to the system Q with small de¬ 
lays. First we want to consider the bifurcations of the 
trivial rest point. From now on, the bifurcation parame¬ 
ters are coupling strength c and time delay t. 


3.2.1. Bifurcations of the trivial rest point 

The type of bifurcations of the rest points determines 
the most important neuro-computational features of the 
neuron. For the system Q neurons are excitable be¬ 
cause the trivial resting state is near a Hopf bifurca¬ 
tion, i.e., near a transition from quiescence to spiking. 
Hence, by increasing c and t, a branch of Hopf bifurca¬ 
tion emanates from c - 0.3974 and t = 0 in (c, t) plane, 
and continues until c - 0.6285 and t = 0.5219 (first br 
of Hopf in Fig. |^, in which the Hopf curve meets the 
curve corresponding to the pitchfork bifurcation of the 
trivial rest point. In this situation the second conditions 
of Lemma 2.5 are satisfied and a Hopf-pitchfork bifur¬ 
cation occurs, see Fig. Since the Hopf bifurcation 
occurs before pitchfork bifurcation, the two nontrivial 
rest points are both unstable. To illustrate what we have 
stated above, we can fix the parameter c and study the 
dynamic changes according to t, see Fig. As it is 


shown in Fig. Ba), for c - 0.5 and r = 0.1 the neu¬ 
rons oscillate almost synchronized, in which the phase 
shift between oscillation of two neurons is close (but not 
equal) to zero. By increasing the parameter t the ampli¬ 
tude of the periodic solutions decrease until the param¬ 
eter T reaches the first branch of Hopf bifurcation. Due 
to the super-critical Hopf bifurcation the stable limit cy¬ 
cle disappears and the trivial rest point becomes stable. 
Fig. Qc). By further increasing the parameter t, the 
neurons start to spike in an almost anti-phase manner, 
in which the phase shift between oscillation of two neu¬ 
rons is close (but not equal) to n. Fig. |5d). We will 
see in the next sections that, by increasing the range 
of parameter t other branches of Hopf bifurcation ap¬ 
pear which are purely delay-driven and have no trace in 
system with instantaneous coupling. We should notice 
that, even for small delays, the parameter r can either 
suppress periodic spiking or induce new periodic spik¬ 
ing, depending on the value of the time delay, see Fig. 
1^ Thus, treatment of neural systems with changing the 
time delay can be delicate and challenging. 

3.2.2. Bifurcations of the non-trivial rest points 

As stated in the previous sections there is a pitch- 
fork bifurcation and two new unstable rest points ap¬ 
pear. We want to study the bifurcations of the nontriv- 
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Figure 5: Bifurcation diagram. Blue line is a branch of Hopf bifurca¬ 
tion for nontrivial rest points for small delays. The red line is eight- 
figure saddle homoclinic branch which emanates at c = 1.0545. The 
dashed black line is the branch of fold of limit cycles, which emanates 
at c = 1.0721. The (red) stared point, shows the pitchfork bifurcation 
when r = 0. The solid green line is the branch of pitchfork bifurca¬ 
tion. 


ial rest points. We will show that the small time de¬ 
lay can induce new period windows with the coupling 
strength increasing. In subsection |3.1| we stated that 
there is an eight-figure saddle homoclinic bifurcation 
near the sub-critical Hopf bifurcation, for the system 
with instantaneous coupling. Therefore, in the system 
with small delayed coupling, in addition to the branch of 
Hopf bifurcation for small delays, there exists a delay- 
induced branch of eight-figure saddle homoclinic bifur¬ 
cation, see Fig. Moreover there is a branch of fold 
of limit cycles that emanates at c ^ 1.0721, close to the 
Hopf and eight-figure saddle homoclinic branches. 

In order to fully understand the bifurcation diagram in 
Fig. 0 we change the coupling parameter c and study 
the dynamic behavior of the system. 

(1) For T = 0.12 and c - 0.7 the bifurcation parameters 
are located below the sub-critical Hopf branch of the 
non-trivial rest point, region A in Fig|^ In this region 
the system is mono-stable and there exists only a stable 
periodic solution. Fig. [^al,a2). 



Figure 6: Left and right figures show time series of vi(f) (solid red 
curve) and V 2 {t) (dashed blue curve) with initial conditions IC = 
[1.3; 1.5; 1.4; 1] and 1C = [0.05; 0.03; 0.04; 0.2] respectively, (al) and 
(a2) Periodic spiking for parameters in region A of Fig. (bl) and 
(b2) Periodic spiking and rest state for parameters in region B of Fig. 

(cl) and (c2) Periodic spiking and rest state for parameters in re¬ 
gion C of Fig. (d 1) and (d2) Rest states for parameters in region D 
of Fig. 1^ 


(2) For T = 0.12 and c = 1, the bifurcation parameters 
are located above the sub-critical Hopf branch of the 
non-trivial rest point and below the homoclinic branch, 
region B in Fig. 0 In this region, the system has a sta¬ 
ble limit cycle and a pair of stable nontrivial rest points. 
With different initial conditions one can reach the stable 
limit cycle and start to spike or tend to the nontrivial rest 
points. Fig. 0^bl,b2). 

(3) For T - 0.12 and c - 1.08 the bifurcation parame¬ 
ters are located above the homoclinic branch and below 
the fold of limit cycles branch, region C in Fig. 0 In 
this region the behavior of the neurons doesn’t change 
and is like what stated above. Fig. 0^cl,c2). 

(4) For T = 0.12 and c = 1.1 the bifurcation parameters 
are located above the fold branch of limit cycles, region 
D in Fig. 0 and only non-trivial rest points are stable 
states in this region, see Fig. 0^dl,d2). 


4. Purely delay driven dynamics 


In this section, in addition to the bifurcations which 
was stated in previous sections, we will study purely 
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Figure 7: Hopf branches of the trivial rest point for a = 0.55, = 

1.128, and b 2 = 0.58. The solid green line is the branch of pitchfork 
bifurcation. 


Figure 8: Real parts of eigenvalues for the trivial rest point vs changes 
of the parameter r. We set parameters a = 0.55,= 1.128, b 2 = 
0.58, and c = 0.2. 


delay driven bifurcations which have no trace in system 
with instantaneous coupling. 


4.1. Bifurcations of the trivial rest point 

We start the study of possible bifurcations of the triv¬ 
ial rest point with finding branches of Hopf bifurcation, 
see Fig. The Hopf branches are the roots of equation 
0. which some of them are shown in Fig. It is not 
hard to see that for parameters c < 0.1016, the condi¬ 


tions of Lemma 2.3 hold and the origin is a stable rest 
point for r > 0. 

As the next step, we increase the parameter c, to meet 
the Hopf branches. For instance for c — 0.2, it is easy to 
check that the conditions of Lemma |2.2[ i) hold, hence 
0 has positive roots. 

We obtain = 0.5739 and = 0.7718, thus. 


0)1 = 0.7575 and 
By 0, 0, and 0 we have. 




5.6096-I-27;r y, 
’’^2 


‘ 2(0.7575) 

j = 0,1,2,... 


a>2 = 0.8785. 


2.869 -I- 2 jn 
2(0.8785) ’ 


(9) 


We should notice that for the above parameters, the con¬ 
ditions of Lemma |23| ii) are satisfied, hence all roots of 
0 have negative real parts, when t e [0, tq), and the 
neuron is at rest, see Fig. @ and Fig. 0a) . By Eq. 
(|^ we have tq = = 1.63. To understand the bi¬ 

furcation diagram in Fig. 0 for c - 0.2 we want to 
change the parameter r and study the dynamic behav¬ 
ior of the system. By equations (01 the first Flopf bi¬ 
furcation occurs for Tf\ see Fig. p and Fig. 8 Ac¬ 
tually, it is a super-critical Hopf bTurcation, thus the 


neuron’s stability switches from resting state to periodic 
spiking, and a small limit cycle branches from the rest 
point. Fig. I^b). Let T 2 = the resulting limit cy¬ 

cle is T-periodic, with T being near T 2 - There¬ 
fore, the period of the small limit cycle attractor ap¬ 
pearing via super-critical Hopf bifurcation is finite and 
positive. Thus, the frequency of oscillations is nonzero, 
but their amplitudes are small. Such periodic oscilla¬ 
tions of small amplitude are not linked to neuronal fir¬ 
ing, but must rather be interpreted as spontaneous sub¬ 
threshold oscillations. The amplitude of the oscillation 
grows with the parameter t up to a limit. Fig. 0c), 
and after that it decreases. Fig. 0d), and when the pa¬ 
rameter reaches t = = 3.7, the super-critical Hopf 

bifurcation occurs and the limit cycle disappears. Fig. 
0e), hence again the trivial rest point becomes stable. 
We should notice that in the stated region the period 
of the periodic solution changes and near the parameter 
T = = 3.7 reaches Ti, which Ti = ^. By increasing 

the parameter r the trivial rest point remains stable for 
T 6 (3.7,5.2). When r = — 5.2, once more again the 

super-critical Hopf bifurcation occurs and the behavior 
of the system changes qualitatively from silence mode 
to 7 ’ 2 -periodic spiking. Fig. 0f). By increasing t the 
above scenario repeats and transitions from silence to 
firing and vice versa occur frequently until = 12.36, 
Fig. |^g,h,k). For parameters t > the trivial rest 
point always has a pair of eigenvalues with positive real 
part, so the trivial rest point remains unstable and the 
neurons spike. In this situation, Hopf bifurcations cause 
the occurrence of a new pair of purely imaginary eigen¬ 
values, Fig. 0 Moreover, there is an interesting ob- 
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Figure 9: Phase portraits of delay-induced oscillations, when c = 0.2. (a) The origin is a stable rest point for r = 1.5 < r®. (b) Small periodic 
oscillation for r = 1.8 £ (t®, t®). (c) Large periodic oscillation for r = 2.5 £ (d) The amplitude of the neuron decreases for r = 3.5. 

(e) The neuron is at rest for t = 4 £ '). (f) Periodic oscillation for t = 5.3 > *. (g) Again, the neuron shows a small periodic oscillation 

for T = 6 £ (Tj \ (h) The neuron shows a large periodic oscillation, (k) The amplitude of the neuron decreases for r = 8.2. 


servation about almost anti-phase and almost synchro¬ 
nized solutions. Actually, almost synchronized and al¬ 
most anti-phase activities of the coupled neurons can be 
achieved in some parameter ranges related to their bifur¬ 
cation transition. If we consider periodic solutions, we 
can see almost anti-phase activities for r e t®), 
in which the period changes from T 2 to Ti, and almost 
synchronized activities for t e 
ilarly the period changes from T 2 to Ti 


, in which sim- 


Fig. 10 For 


an explanation we should consider the eigenspace of 
the linearized system. The eigenfunction coiTesponding 
to the eigenvalue ia )2 is M 2 , M 3 , M 4 ), 

where ui - U 2 {b\ -t- i(x> 2 ), M 3 = U 2 e'‘^^'^'^^ \~(a-ioj 2 )(bi+ia> 2 ) ^ 

and M 4 = It is easy to check that 


when i is an even number (including zero), - 

cos(y -h jn) - 1 - !sin(t^ - 1 - jn) lies in the first quadrant, 
so if M 3 and M 4 lie in the first quadrant, mi and M 2 will 
lie in the fourth quadrant and vice versa. Therefore, two 
neurons have almost anti-phase oscillations. When j is 


an odd number, = cos(^ + jn) + !sin(^ + jn) 

lies in the forth quadrant. Thus, mi, M 2 , M 3 , and M 4 lie in 
the same quadrant. Therefore, two neurons have almost 
synchronized oscillations. As a result, it is seen that 
time delay can influence dynamical behavior of coupled 


neurons, such as changing between altering almost syn¬ 
chronized or almost anti-phase oscillations. The stated 
phenomenon is a consequence of double Hopf bifurca¬ 
tion which is of codimension-two ll58l . 

4.2. Double Hopf bifurcation 

In the previous sections we reviewed Hopf bifurca¬ 
tions of the trivial rest point and the multiple stability 
regions which are formed in the (t, c) plane, and im¬ 
ply multi-time switches between the stable and unstable 
states of the zero solution. In this section we want to 
study other possible bifurcations of the trivial rest point, 
Fig.[n] There are some points in the bifurcation dia¬ 
gram of the trivial rest point. Fig. [m which can be ob¬ 
tained from the intersection of the two branches of Hopf 
bifurcation. It is possible that the characteristic equation 
at a rest point has eigenvalues with strict nonzero real 
parts except two pairs of purely imaginary eigenvalues; 
in (t, c) plane it happens when two branches of Hopf 
bifurcation denoted by ti and T 2 correspond to wj and 
C 02 respectively, cross each other in a point (T^/,,Crf/,). 
This situation is called a double-Hopf bifurcation ll58l . 
For constants let +/aJi and ±iu 2 be the two 

pairs of purely imaginary eigenvalues of the characteris¬ 
tic equation and all other eigenvalues have non zero real 
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Figure 10: Delay-induced oscillations and different modes of oscillation for solutions depicted in Fig. (a) The origin is a stable rest point for 
r = 1.5 < (b) The neurons oscillate in an anti-phase manner when r = 1.8 6 (c) The neurons oscillate in an almost anti-phase 

manner, when r = 2.5 6 (t 2 ^\ (d) Two coupled neurons oscillate in an almost anti-phase manner, r = 3.5. (e) The neurons are at rest for 

T = 4 e (f) The neuron shows periodic oscillation for r = 5.3 > T 2 ^ (g) Again, the neurons show a small periodic oscillation but 

synchronous for r = 6 6 (h) The neurons show large periodic oscillations which are almost synchronized, (k) The amplitude of the 

neurons decrease for r = 8.2. 



Figure 11: Branches of Hopf (blue solid lines), torus (red dotted lines), and pitchfork cycle (black dashed lines) bifurcations of the trivial rest point. 
The solid green line is the branch of pitchfork bifurcation. 
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Figure 12: The different regions ai'ound a double-Hopf bifurcation 
of the non-trivial rest point which is related to second and third Hopf 
branches. 


Figure 13: Two different limit cycles with different frequencies which 
exist when the parameters are in region D. One of them is correspon¬ 
dent to synchronized solutions and another one is con'espondent to 
anti-phase solutions, Fig. El 


parts. The existence of double-Hopf bifurcation leads 
to different modes of spiking in different regions around 
the double-Hopf point, regions C,D, and E in Fig. [T^ 

It follows from the previous paragraph that different val¬ 
ues of delay and coupling parameters of the system Q. 
can induce two different values of frequencies corre¬ 
sponding to the second and third Hopf bifurcation. The 
existence of a double-Hopf bifurcation can be detected 
by the local stability analysis. In this section, we rep¬ 
resent the typical dynamical behaviors and the bifurca¬ 
tion sets in the neighborhood of one of the double-Hopf 


bifurcation points of Fig. 11 First, time delay and cou¬ 
pling strength are chosen from the region A of Fig. 12 
In this region the trivial rest point is locally asymptot¬ 
ically stable. By increasing the parameter c, the trivial 
rest point loses its stability in result of the super-critical 
Hopf bifurcation, where the parameters encounter the 
second Hopf bifurcation branch, and enter region B of 
Fig. [T^ The system periodically oscillates with the u)i 
frequency. On the other hand, when the parameter t is 


increased from region A to F of Fig. 12 the parameters 
will encounter the third Hopf bifurcation branch related 
to 0)2 frequency. This implies that the trivial rest point 
loses its stability through a super-critical Hopf bifurca¬ 
tion and the system exhibits the oscillatory solution with 
0)2 frequency. If we start from region B with a fixed 
parameter c and increase the parameter t, the parame¬ 
ters encounter the third super-critical Hopf bifurcation 
branch and enter region C. Since the trivial rest point is 
unstable, an unstable limit cycle appears. In this situ¬ 
ation, the system is mono-stable and exhibits the oscil¬ 
latory solution with o)\ frequency. By further increas¬ 


ing the the parameter t from region C to D, the unsta¬ 
ble periodic solution with wa frequency becomes stable 
through the sub-critical torus bifurcation. It suggests 
that the system behavior changes to co-existence of two 
different periodic solutions with different frequencies. 
Actually for parameters in region D there are two dif¬ 


ferent stable limit cycles as depicted in Fig. 13 and the 


system is bi-stable. We should notice that in this region 
the behavior of the system can be different when we start 
with different initial conditions and two neurons with¬ 
out changing the parameters can show oscillatory solu¬ 
tions with two different frequencies, and different types 
of synchronization. Fig. 14 We should emphasis that 


there is an unstable torus in region C. Torus bifurcation 
occurs when a complex-conjugate pair of Floquet mul¬ 
tipliers crosses the unit circle Il57l . In the above stated 
torus, a complex-conjugate pair of Floquet multipliers 
moves from the outside to the inside of the unit circle, 
Fig[^ When the time delay is increased from region 
D to E of Fig. [T^ the stable periodic solution which ap¬ 
peared through the second branch of Hopf bifurcation 
loses its stability through the sub-critical torus bifurca¬ 
tion. Hence, in region D the system becomes mono¬ 
stable. For the stated torus a complex-conjugate pair of 
Floquet multipliers moves from inside to the outside of 
the unit circle. Fig. 

As we showed in the figures, Hopf bifurcation gives rise 
to the periodical activity, while double-Hopf bifurcation 
results in two types of oscillations with different fre¬ 
quencies and different synchrony, and a small change 
in the time delay can result in different patterns of activ¬ 
ity. This is particularly useful when using the time delay 
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Figure 14: Different modes of spiking due to different initial condi¬ 
tions, here c - 0.325 and r = 4.7756. (a) Almost anti-phase oscil¬ 
lations when IC=(0.1,0.3,0.4,0.2). (b) Almost synchronized oscilla¬ 
tions when IC=(1,0.9,0.8,0.7). 



Figure 17: Bifurcations of the non-trivial rest points. The solid blue 
lines are Hopf bifurcation branches. The dash-dot black line is the 
pitchfork cycle bifurcation. The dashed black line is fold bifurcation 
of limit cycles. The dotted red lines are torus bifurcation branches. 
The solid red line is fold of limit cycles branch. 



Figure 15: (Left) Floquet multipliers before first torus bifurcation. 
(Right) Floquet multipliers after first torus bifurcation. 



Figure 16: (Left) Floquet multipliers before second torus bifurcation. 
(Right) Floquet multipliers after second torus bifurcation. 


to control bifurcations and oscillatory solutions. 

As it is depicted in Fig. 0 there are some branches of 
pitchfork cycle bifurcation. We will explain their roles 
on the system’s dynamic in the next section. 


4.3. Bifurcations of the non-trivial rest points 


In section 4.1 according to the bifurcations of the 


trivial rest point, we showed different stability regions in 
the (t, c) plane, which implies multi-time switches be¬ 
tween the stable and unstable states of trivial solution. 
Similarly we can see such multiple stability regions ac¬ 
cording to bifurcations of non-trivial rest points. In this 
way, we will study the possible bifurcations of the non¬ 
trivial rest points too, see Fig. 17 In the following. 


we restrict ourselves to a fixed value of c and study the 
delay induced changes of the dynamics according to all 
possible bifurcations of the system. Since in the study 
of the dynamics of the system ([T]i, the bifurcations of 
the origin and non-trivial rest points are both indispens¬ 
able, we draw entire bifurcation diagram which were 


depicted in Fig. 17 and Fig. 11 in a single figure, see 
Fig. For fixed parameter c - 1.03, we change the 
parameter t, and investigate the behavior of the system. 
When T = 0.05, the parameters are located in region 1 
of Fig. [T^and the dynamics is the same as in region B 
in Fig. ppTherefore, the system has a stable limit cycle 
related to synchronous oscillations, and a pair of stable 
non-trivial rest points. By increasing the parameter r. 
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Figure 18: Total bifurcation diagram of the system Blue solid lines are Hopf branches of non-trivial rest points. Magenta dash-dot lines are 
Hopf branches of trivial rest point. Black dashed lines are pitchfork bifurcations of limit cycles. Red dotted lines are toms bifurcation branches. 
The other branches are introduced in Fig. 


for example r = 0.5, the parameters are located in re- 


18 By further increasing of the parameter r, the third 


gion 2 of Fig. 18 and the dynamics is as region D in 


Fig. 1^ Since the parameters are located after the branch 
of fold of limit cycles, only non-trivial rest points are 
stable. By a little change of the parameter r the second 
branch of super-critical Hopf bifurcation of the trivial 
rest point occurs, and the trivial rest point which was 
unstable gives birth to a small-amplitude unstable limit 


cycle. Therefore, in region 3 of Fig. 18 the two non¬ 
trivial rest points remain stable. By increasing the pa¬ 
rameter T a sub-critical pitchfork bifurcation of limit cy¬ 
cles occurs, the stated unstable limit cycle changes sta¬ 
bility, becomes stable and produces two unstable limit 
cycles. In result, in this region, according to the initial 
conditions the neurons spike in an anti-phase manner, or 
remain in the rest states (non-trivial rest points), region 
4 in Fig. 18 By more changes of r the delay-induced 


limit cycle stay stable and never disappear, but each of 
the unstable limit cycles goes towards a non-trivial rest 
point and their amplitude decrease. When t reaches the 
second branch of sub-critical Hopf bifurcation of non¬ 
trivial rest points in Fig. 17 the small unstable limit 


cycles shrinks to stable non-trivial rest points and make 
them lose stability. Hence, the system becomes mono¬ 
stable and the neurons start to spike, region 5 in Fig. 


branch of super-critical Hopf bifurcation for trivial rest 
point appears. Fig. El and the trivial rest point, which 
is unstable of Morse index 3, gives birth to a small- 
amplitude unstable limit cycle with index 3. Therefore, 


in region 6 of Fig. 18 the system remains mono-stable. 
By increasing the parameter t again a sub-critical pitch- 
fork bifurcation of limit cycles occurs, and the index-3 
unstable limit cycle changes to index-2 unstable limit 
cycle. Hence, in region 7 of Fig. 18 the system re¬ 


mains mono-stable. If we increase the parameter r, a 
sub-critical torus bifurcation occurs, and the index-2 un¬ 
stable limit cycle changes stability and becomes stable. 
The dynamics in region 8 of Fig. is the same as dy¬ 
namics in region D of Fig. Therefore, the system 
behavior changes to co-existence of two different peri¬ 
odic solutions with different frequencies and different 
synchrony. Actually, for parameters in region 8 there 
are two different stable limit cycles, one of them is re¬ 
lated to almost synchronous oscillations and the other 
one is related to almost anti-phase activities. The dy¬ 
namics in region 9 of Fig[^ is really interesting and 
complicated. By changing parameters from region 8 
to 9, a sub-critical Hopf bifurcation for non-trivial rest 
points occurs and makes the unstable non-trivial rest 
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Figure 19: Different dynamical activities in region 9 of Fig. |18| Whit 
different initial conditions the system reaches different stable states. 
Flere c = 1.08 and r = 3.9. 


points stable. Hence, in region 9, there are two dif¬ 
ferent periodic solutions with different frequencies, and 
also a pair of stable rest points, see Fig. Related to 
double-Hopf bifurcation point which is stated below the 
region 9, there are branches of super-critical torus bifur¬ 
cations. Since the Hopf bifurcations of the non-trivial 
rest points are sub-critical, the limit cycles which are 
related to these bifurcations and also the created torus 
are unstable. If we increase t from region 9 to 10 again 
a sub-critical Hopf bifurcation occurs for non-trivial rest 
points and makes them unstable. Hence in region 10 of 
Fig. 18 there are two different stable limit cycles, and 
rest points are unstable. 


4.4. Strong coupling 

For the system ([T]l with instantaneous coupling (r = 
0 ) the characteristic equation at the trivial rest point is 

P{c, t) = / H- AA^ + BA^ + CA + D ^0, (10) 

where A = b\+b2—2a, B = b\b2—2a{b\+b2)+a^—c^+2, 
C - (a^ + l)(bi + ^ 2 ) - 2abib2 — 2a — c^{bi + ^ 2 ), and 
D - (a^- c^)bib 2 - a(bi + ^ 2 ) + 1. According to Lemma 


least one positive root. It is easy to check that for chosen 
values of the parameter in our paper, when c > 0.6265, 
D is negative and the characteristic equation has at least 
a positive root. Indeed, for c > 0.6265 the trivial rest 
point is unstable of index 1. Moreover, in the system 
with delayed coupling ([T]), when pitchfork bifurcation 
occurs the trivial rest point becomes unstable of index 
k, where k is an odd number. Since the delay driven 
bifurcations such as Hopf bifurcation change the rest 


2.1 ifD < 0, then the characteristic equation (|10[) has at 


Figure 20: Different limit cycles with different frequencies which exist 
in region X of Fig. |18| The dash-dot black line is correspondent to 
synchronized solutions and the others are correspondent to anti-phase 
solutions. 



Figure 21: The dash-dot black line is correspondent to synchronized 
solutions, the solid blue line is correspondent to anti-phase solutions 
with highter frequency, and the dotted red line is correspondent to 
anti-phase solutions with lower frequency. 


point of index k to a rest point of index k ±2, the in¬ 
dex of the rest point remains odd, and the trivial rest 
point remains unstable. Although the trivial rest point is 
unstable for all the parameter space in Fig. the influ¬ 
ence of it’s Hopf bifurcations on the dynamics isn’t neg¬ 
ligible. Our interest is to determine if there are some 
stable solutions, and synchrony properties of oscillatory 
activities of the neurons. In this way, for hxed parame¬ 
ter c = 1.2, we change the parameter t, and investigate 
the behavior of the system. When the parameters are 
located in region I of Fig. 18 the dynamics is the same 
as region 2 of Fig. 18 Therefore, only the two non¬ 
trivial rest points are stable and they remain stable for 
all values of parameter r in parameter space of Fig. 

By further change of the parameter t the second super- 
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critical Hopf bifurcation of the trivial rest point occurs, 
and the trivial rest point which was unstable of index 1 
gives birth to a small-amplitude unstable limit cycle of 
index 1. Therefore in region II of Fig. [T^ the two non¬ 
trivial rest points are stable, and there exists an unstable 
limit cycle. By a little change of the parameter t a sub- 
critical pitchfork bifurcation of limit cycles occurs, the 
stated unstable limit cycle changes stability, becomes 
stable and produces two unstable limit cycles. In result, 
in region III of Fig. according to the initial condition 
the neurons spike in an almost anti-phase manner, or re¬ 
main in the rest state (non-trivial rest points). By further 
increasing of the parameter t, the third branch of super¬ 
critical Hopf bifurcation of the trivial rest point appears. 
Fig. El and the trivial rest point, which is unstable of 
index 3, gives birth to a small-amplitude unstable limit 
cycle with index 3. Therefore, the activity of the neu¬ 
rons in region IV of Fig. 18 remains similar to those in 
region III. By increasing the parameter t, again a sub- 
critical pitchfork bifurcation of limit cycles occurs, and 
the index-3 unstable limit cycle changes to index-2 un¬ 
stable limit cycle. Hence, in region V of Fig. 18 the 
activity of the neurons is the same as region IV. If we 
increase the parameter t, a sub-critical torus bifurcation 
occurs, and the index-2 unstable limit cycle changes sta¬ 
bility and becomes stable. This stable solution is related 
to almost synchronized oscillations. The dynamics in 
region VI of Fig. 18 is the same as dynamics in region 
9, hence there are two different periodic solutions with 
different frequencies and different synchrony, and also a 
pair of stable rest points, see Fig. As t is increased 
from region VI to VII, another branch of super-critical 
Hopf bifurcation of the trivial rest point appears, and the 
trivial rest point which is unstable of index 5 becomes 
unstable of index 7. This Hopf bifurcation results in the 
creation of a small-amplimde unstable limit cycle with 
index 5. By a little change of the parameter t, a sub- 
critical pitchfork bifurcation of limit cycles occurs, and 
the index-5 unstable limit cycle changes to index-4 un¬ 
stable limit cycle, region VIII in Fig. [T^ By further 
change of the parameter t, a torus bifurcation occurs 
and the index-4 unstable limit cycle changes to index-2 
unstable limit cycle, region IX in Fig. The activity 
of the neurons is not affected qualitatively by this bifur¬ 
cation. When the parameter r is increased from region 
IX to X, the index-2 unstable limit cycle changes stabil¬ 
ity and becomes stable through a torus bifurcation. The 
global dynamics in region X is characterized by three 
different stable limit cycles which are depicted in Fig. 
[20] These oscillatory solutions consist of two almost 
anti-phase and one almost synchronized solutions. Two 
almost anti-phase solutions are recognizable throughout 


their different frequencies, as is depicted in Fig. |2T| 
Therefore, the system as a whole is multi-stable, two 
non-trivial rest points and three limit cycles. As the pa¬ 
rameter T changes from region X to XI, the limit cycle 
which had become stable in region III and is of lower 
frequency, changes stability and becomes unstable of 
index 2. Hence, in this region the dynamics could corre¬ 
spond to either synchronized and anti-phase spiking, or 
quiescence. We should remind that, for narrow regions 
of delay parameter, by the mechanism stated above and 
in the previous section, the parameter t can either sup¬ 
press periodic spiking, or induce new ones with different 
frequencies and different synchrony. 

5. conclusion 

As a typical representation of excitable systems, FHN 
model is widely used to research some nonlinear phe¬ 
nomena in the nervous system. In this paper, possible 
bifurcations were described in two synaptically coupled 
FHN neurons with the time delay. We showed that with 
the coupling strength changing, two coupled nonidenti¬ 
cal neurons can exhibit rich bifurcation behaviors. We 
also showed that the dynamics in two coupled FHN neu¬ 
rons is drastically changed due to the effect of time de¬ 
lay. We used theory of dynamical systems to describe 
typical behaviors, such as quiescence, almost synchro¬ 
nized and anti-phase periodic solutions, or coexistence 
of these states. It illustrates that the time delay can either 
suppress periodic activity or induce new ones, depend¬ 
ing on the value of the time delay and kind of bifurca¬ 
tion. All possible Hopf, double Hopf, pitchfork, pitch- 
fork of limit cycles, homoclinic, fold of limit cycles, 
and torus bifurcations of our system are derived, and to¬ 
tal bifurcation diagrams for trivial and non-trivial rest 
points are reported. Actually this kind of bifurcation 
study is a serious task for the prediction and detection 
of all possible dynamical behaviors. Missing some kind 
of solutions is possible, if we only aid simulation study. 
We should emphasis that our investigation provides a 
detailed scenario for the Hopf bifurcation diagrams ap¬ 
peared in similar studies. These results are interesting 
from the point of view of applications, since our generic 
model is representative for a wide range of real-world 
systems. In the future research we will classify possi¬ 
ble activities for two coupled nonidentical FHN neurons 
without inherent symmetry. 
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